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3 , Abstract 

CZ i 

In most books the Delaunay and Lagrange equations for the orbital elements are de- 
rived by the Hamilton-Jacobi method: one begins with the two-body Hamilton equations 
in spherical coordinates, performs a canonical transformation to the orbital elements, and 
obtains the Delaunay system. A standard trick is then used to generalise the approach 
to the N-body case. We re-examine this step and demonstrate that it contains an im- 
plicit condition which restricts the dynamics to a 9(N-l)-dimensional submanifold of the 
12(N-l)-dimensional space spanned by the elements and their time derivatives. The tacit 
condition is equivalent to the constraint that Lagrange imposed "by hand" to remove the 
excessive freedom, when he was deriving his system of equations by variation of parame- 
ters. It is the condition of the orbital elements being osculating, i.e., of the instantaneous 
ellipse (or hyperbola) being always tangential to the physical velocity. Imposure of any 
supplementary condition different from the Lagrange constraint (but compatible with the 
equations of motion) is legitimate and will not alter the physical trajectory or velocity 
(though will alter the mathematical form of the planetary equations). 
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This freedom of nomination of the supplementary constraint reveals a gauge-type 
internal symmetry instilled into the equations of celestial mechanics. Existence of this 
internal symmetry has consequences for the stability of numerical integrators. Another 
important aspect of this freedom is that any gauge different from that of Lagrange makes 
the Delaunay system non-canonical. In a more general setting, when the disturbance de- 
pends not only upon positions but also upon velocities, there is a "generalised Lagrange 
gauge" wherein the Delaunay system is symplectic. This special gauge renders orbital 
elements that are osculating in the phase space. It coincides with the regular Lagrange 
gauge when the perturbation is velocity-independent. 



1 Euler and Lagrange 
1.1 The history 

The planetary equations, which describe the evolution of the orbital elements, constitute the 
cornerstone of the celestial mechanics. Description of orbits in the language of Keplerian 
elements (rather than in terms of the Cartesian coordinates) is not only physically illustrative 
but also provides the sole means for analysis of resonant interactions. These equations exist 
in a variety of equivalent forms (those of Lagrange, Delaunay, Gauss, Poincare) and can be 
derived by several different methods. 

The earliest sketch of the method dates back to Euler's paper of 1748, which addresses 
the perturbations exerted upon one another by Saturn and Jupiter. In the publication on the 
Lunar motion, dated by 1753, Euler derived the equations for the longitude of the node, Q , 
the inclination, % , and the quantity p = a (1 — e 2 ) . Time derivatives of these three elements 
were expressed through the components of the disturbing force. Sixty years later the method 
was amended by Gauss who wrote down similar equations for the other three elements and, 
thus, obtained what we now call the Gauss system of planetary equations. (The history of 
this scientific endeavour was studied by Subbotin (1958), who insists that the Gauss system 
of planetary equations should rather be called Euler system.). A modern but still elementary 
derivation of this system belongs to Burns (1976). 

In his memoires of 1778, which received an honourable prize from the Academie des Sciences 
of Paris, Lagrange employed the method of variation of parameters (VOP) to express the time 
derivatives of the orbital elements through the disturbing functions' partial derivatives with 
respect to the Cartesian coordinates. In his memo ire of 1783 Lagrange furthered this approach, 
while in Lagrange (1808, 1809, 1810) these equations acquired their final, closed, shape: they 
expressed the orbital elements' evolution in terms of the disturbing potentials' derivatives 
with respect to the elements. Lagrange's derivation rested upon an explicit imposure of the 
osculation condition, i.e., of a supplementary vector equation (called the Lagrange constraint) 
which guaranteed that the instantaneous ellipses (in the case of bound motions) or hyperbolae 
(in the case of flyby ones) are always tangential to the physical trajectory. Though it has 
been long known (and, very possibly, appreciated by Lagrange himself) that the choice of the 
supplementary conditions is essentially arbitrary, and though a couple of particular examples 
of use of non-osculating elements appeared in the literature (Goldreich 1965; Brumberg et al 
1971; Borderies & Longaretti 1987), a comprehensive study of the associated freedom has not 
appeared until very recently (Efroimsky 2002, 2003). 

In the middle of the 19th century Jacobi applied a canonical-transformation-based procedure 
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(presently known as the Hamilton- Jacobi approach) to the orbital dynamics, and offered a 
method of deriving the Lagrange system. This technique is currently regarded standard and is 
offered in many books. Though the mathematical correctness of the Hamilton- Jacobi method 
is beyond doubt, its application to celestial mechanics contains an aspect that has long been 
overlooked (at least, in the astronomical literature). This overlooked question is as follows: 
where in the Hamilton- Jacobi derivation of the planetary equations is the Lagrange constraint 
tacitly imposed, and what happens if we impose a different constraint? This issue will be 
addressed in our paper. 



1.2 The gauge freedom 

Mathematically, we shall concentrate on the N-body problem of celestial mechanics, problem 
that for each body can be set as 

? + 4 - = AF , (1) 

AF being the disturbing force that vanishes in the (reduced) two-body case, r being the 
position relative to the primary, and /i standing for G(m p i ane t + m sun . A solution to the 
unperturbed problem is a Keplerian ellipse (or hyperbola) 

r = i(C u ...,C 6 ,t) (2) 

parametrised by six constants (which may be, for example, the Kepler or Delaunay elements). 
In the framework of the VOP approach it gives birth to the ansatz 

r = F^t),...,^*),*) , (3) 

the " constants" now being time-dependent and the functional form of f remaining the same as 
in (J2J). Substitution of (j3J) into (JTJ) results in three scalar equations for six independent functions 
Ci(t) . In order to make the problem defined, Lagrange applied three extra conditions 

V — — = (4) 

V dQ dt ' 1 } 

that are often referred to as "the Lagrange constraint." This constraint guarantees osculation, 
i.e., that the functional dependence of the perturbed velocity upon the "constants" is the same 
as that of the unperturbed one. This happens because the physical velocity is 

where g stands for the unperturbed velocity that emerged in the two-body setting. This 
velocity is, by definition, a partial derivative of f with respect to the last variable: 

g((7i, C 6 , t) = — . (6) 

This choice of supplementary conditions is convenient, but not at all necessary. A choice of 
any other three scalar relations (consistent with one another and with the equations of motion) 
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will give the same physical trajectory, even though the appropriate solution for non-osculating 
variables C{ will differ from the solution for osculating ones. 

Efroimsky (2002, 2003) suggested to relax the Lagrange condition and to consider 

$ being an arbitrary function of time, "constants" C{ and their time derivatives of all orders. 
For practical reasons it is convenient to restrict $ to a class of functions that depend upon the 
time and the "constants" only. (The dependence upon derivatives would yield higher-than- 
first-order time derivatives of the in the subsequent developments, which would require 
additional initial conditions, beyond those on r and r , to be specified in order to close the 
system.) Different choices of $ entail different forms of equations for Cj and, therefore, 
different mathematical solutions in terms of these "constants." A transition from one such 
solution to another will, though, be a mere reparametrisation of the orbit. The physical orbit 
itself will remain invariant. Such invariance of the physical content of a theory under its 
mathematical reparametrisations is called gauge symmetry. On the one hand, it is in a close 
analogy with the gradient invariance of the Maxwell electrodynamics and other field theories. 
On the other hand, it illustrates some general mathematical structure emerging in the ODE 
theory (Newman & Efroimsky 2003). 

If the Lagrange gauge (HJ) is fixed, the parameters obey the equation 

[C n Cj] standing for the unperturbed (i.e., defined as in the two-body case) Lagrange brackets: 

[Cn Cj] = dC n dCj ~ dC~ dG n ' (9) 

To arrive at formula (JSJ), one should, according to Lagrange (1778, 1783, 1808, 1809), differenti- 
ate (J3J), insert the outcome into (QJ, and then combine the result with the Lagrange constraint 
(J1J. (In the modern literature, this derivation can be found, for example, in Brouwer and 
Clemence (1961), Efroimsky (2002, 2003), Newman and Efroimsky (2003), Efroimsky and 
Goldreich (2003).) 

In the simplest case the perturbing force depends only upon positions and is conservative: 
AF = dR(r)/dr . Then the right-hand side of (JHJ) will reduce to the partial derivative 
of the disturbing function R(r) with respect to C n , whereafter inversion of the Lagrange- 
bracket matrix will entail the Lagrange system of planetary equations (for Cj being the Kepler 
elements) or the Delaunay system (for the parameters chosen as the Delaunay elements). 

As explained in Efroimsky (2003), in an arbitrary gauge $ equation (JSJ) will generalise to 
its gauge-invariant form 

L [Cn C,] w - ^- AF - _ — - _ $ , (10) 

the Lagrange brackets [C n Cj] being still defined through (|9~|). If we agree that $ is a 
function of both time and the parameters Ci , but not of their derivatives, then the right-hand 
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side of fllOj) will implicitly contain the first time derivatives of Ci . It will then be reasonable 
to move them into the left-hand side. Hence, ()10|) will be reshaped into 



E ( [CnCj] + da d ^ ) dt 3 da AF m ^ $ (11) 



This is the general form of the gauge-invariant perturbation equations of celestial mechanics, 
that follows from the VOP method, for an arbitrary disturbing force AF(r, r, t) and under 
the simplifying assumption that the arbitrary gauge function $ is chosen to depend on the 
time and the parameters Cj , but not on their derivatives. 

For performing further algebraic developments of (fTUj) and (JUJ), let us decide what sort of 
interactions will fall within the realm of our interest. On general grounds, it would be desirable 
to deal with forces that permit description in the language of Lagrangians and Hamiltonians. 



2 Delaunay 

2.1 Perturbations of Lagrangians and Hamiltonians 

Contributions to the disturbing force AF generally consist of two types, physical and inertial. 
Inputs can depend upon velocity as well as upon positions. As motivation for this general- 
ization we consider two practical examples. One is the problem of orbital motion around a 
precessing planet: the orbital elements are defined in the precessing frame, while the velocity- 
dependent fictitious forces play the role of the perturbation (Goldreich 1965, Brumberg et al 
1971, Efroimsky & Goldreich 2003). Another example is the relativistic two-body problem 
where the relativistic correction to the force is a function of both velocity and position, as 
explained, for example, in Brumberg (1992) and Klioner & Kopeikin (1994). (It turns out that 
in relativistic dynamics even the two-body problem is disturbed, the relativistic correction act- 
ing as disturbance. This yields the gauge symmetry that will cause ambiguity in defining the 
orbital elements of a binary.) Finally, we shall permit the disturbances to bear an explicit time 
dependence. Such a level of generality will enable us to employ our formalism in noninertial 
coordinate systems. 

^ 2 

Let the unperturbed Lagrangian be r /2 — U(r). The disturbed motion will be described 
by the new, perturbed, Lagrangian 

2 

£ = r — - U(r) + A£(r, ?, t) , (12) 

and the appropriately perturbed canonical momentum and Hamiltonian: 

u dAC . , 

P = ? + — - , (13) 
or 

-> 2 i / o A f \ 2 

H = p?-C = ^- + U + AH , AH = — AC - - ( — ^ . (14) 

The Euler-Lagrange equations written for the perturbed Lagrangian ()12J) are: 

dU 

? = - — + AF , (15) 
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where the disturbing force is given by 



AF 



dAC 
dr 



d_ ( dAC \ 



(16) 



We see that in the absence of velocity dependence the perturbation of the Lagrangian plays the 
role of disturbing function. Generally, though, the disturbing force is not equal to the gradient 
of AC , but has an extra term called into being by the velocity dependence. 

As we already mentioned, this setup is sufficiently generic. For example, it is convenient 
for description of a satellite orbiting a wobbling planet: the inertial forces, that emerge in the 
planet-related noninertial frame, will nicely fit in the above formalism. 

It is worth emphasising once again that, in the case of velocity-dependent disturbances, 
the disturbing force is equal neither to the gradient of the Lagrangian's perturbation nor to 
the gradient of negative Hamiltonian's perturbation. This is an important thing to remember 
when comparing results obtained by different techniques. For example, in Goldreich (1965) 
the word " disturbing function" was used for the negative perturbation of the Hamiltonian. For 
this reason, the gradient of so defined disturbing function was not equal to the disturbing force. 
A comprehensive comparison of the currently developed theory with that offered in Goldreich 
(1965) will be presented in a separate publication (Efroimsky & Goldreich 2003), where we 
shall demonstrate that the method used there was equivalent to fixing a special gauge (one 
described below in subsection 2.3 of this paper). 



2.2 Gauge- invariant planetary equations 

Insertion of the generic force (j!6)l into ()10|) will bring us: 



E Pn Ci 
3 



df dAC 



dt dC n df dC n dt 
If we recall that, for a velocity-dependent disturbance, 

dAC dAC df dAC d'r dAC df 

+ 



df d (-* dAC\ dg 
dr 7 dC n 



dAC d{§ + $) 



dC n df dC n df dC n 
then equality (fT7|) will look like this: 

^ ^ dCj dAC dAC d® 

E[^]— = ^ 



df dC n 



df 



df d 



dC n 



df dC n 



dCr, dt 



d_AC\ 
dr J 



dC n 



dg 

dC„ 



(17) 



(18) 



After subsequent addition and subtraction of (l/2)d(( y d( y AC)/dr) 2 )/dC n in the ri 

— *■ • 

side, the gauge function $ will everywhere appear in the company of + d(AC)/dr 




E [G. Ci] + 



df d (dAC 



dC n dCj \ dr 




(20) 



d 



da, 



. „ 1 (dACY 



dg 

dC n 



df d 



dAC d 



dC ri dt 




dAC' 
dr , 
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the sum in square brackets being equal to — AH. While (|TT|) expressed the VOP method 
in the most generic form it can have in terms of disturbing forces AF(r, r, t), equation ()20j) 
furnishes the most general form in terms of the Lagrangian perturbation AC(r, r, t) (under 

— * 

the simplifying assumption that the arbitrary gauge function <E» is set to depend only upon 
the time and the parameters Cj , but not upon their derivatives). 

The Lagrange brackets in (|19|) are gauge- invariant; they contain only functions f and g 
that were defined in the unperturbed, two-body, setting. This enables us to exploit the well- 
known expressions for the inverse of this matrix. These look most simple (and are either zero 
or unity) in the case when one chooses as the "constants" the Delaunay set of orbital variables. 
As well known, this simplicity of the Lagrange and their inverse, Poisson, brackets of the 
Delaunay elements is the proof of these elements' canonicity in the unperturbed, two-body, 
problem. When only a position-dependent disturbing function R(f) = AC(r) is "turned 
on," the Delaunay elements still remain canonical, provided the Lagrange gauge is imposed. 
This happens because, as is well known (Brouwer & Clemence 1961), the equations of motion 
together with the Lagrange constraint yield, in that case, the following equation: 

HC B AC 

E[CnC,]^= — , AC = AC ...,C 6 ,t)) = R(f(C h ...,C 6 ,t)) ,(21) 

j n 

which, in its turn, results in the standard Delaunay system. 

In our case, though, the perturbation depends also upon velocities; beside this, the gauge 
$ is set arbitrary. Then our equation fl2()(l will entail the gauge-invariant Lagrange-type 
and Delaunay- type systems of equations that are presented in Appendix 1. Interestingly, the 
gauge-invariant Delaunay-type system is, generally, non-symplectic. It regains the canonical 
form only in one special gauge considered below (a gauge which coincides with the Lagrange 
gauge when the perturbation bears no velocity dependence). This can be proven either by 
a direct substitution of that special gauge condition into the gauge-invariant Delaunay-type 
system given in Appendix 1. An easier way would be to fix the gauge already in ([20j). and this 
is what we shall do in the next subsection. 



2.3 The generalised Lagrange gauge: 

gauge wherein the Delaunay-type system becomes canonical 

We transformed (|17|) into (j2~Uj) for two reasons: to single out the negative perturbation of the 
Hamiltonian and to reveal the advantages of the gauge 

d AL , , 

* = - — — , (22) 

OT 

— *■ 

which reduces to $ = for velocity-independent perturbations. The first remarkable peculiar- 
ity of (|22j) is that in this gauge the canonical momentum p is equal to g (as can be seen from 
© and Pjl): 

u - u dAC _ , , 

g = r- $ = r + — ^- = p . (23) 

dr 

We see that in this gauge not the velocity but the momentum in the disturbed setting is the 
same function of time and "constants" as it used to be in the unperturbed, two-body, case. 
Stated differently, the instantaneous ellipses (hyperbolae) defined in this gauge will osculate 
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the orbit in the phase space. For this reason our special gauge (|22|) will be called the 
"generalised Lagrange gauge." 

Another good feature of (J22J) is that in this gauge equation (J2UJ) acquires an especially 
simple form: 



whose advantage lies not only in its brevity, but also in the invertibility of the matrix emerg- 
ing in its left-hand side. As already mentioned above, the gauge invariance of definition © 
enables us to employ the standard (Lagrange-gauge) expressions for [C n Cj\~ l and, thus, to 
get the planetary equations by inverting matrix [C n Cj] in ({Hj). The resulting gauge- invariant 
Lagrange- and Delaunay-type systems are presented in Appendix 1. 

In the special gauge (J22J), however, the situation is much better. Comparing (}2*T]) with (JUJ), 
we see that in the general case of an arbitrary R = AC(r, r, t) one arrives from (J24j) to the 
same equations as from (|21|) . except that now they will contain — AH. instead of R = AC 
These will be the Delaunay-type equation in the generalized Lagrange gauge: 



(25) 



(26) 



(27) 



dL dAH d( - M ) dAH 

~dt ~ d{ - M ) ' dt ~ dL ' 

dG dAH d(-oo) dAH 

~dt ~ d(-u) ' dt ~ dG ' 

dH dAH d(-Sl) dAH 

~dT ~ d(-Q) ' dt ~ dH 

where 

L = /i 1 / 2 a 1 ' 2 , G = v l/2 a l l 2 (l - e 2 )''" , H = fi 1/2 a 1 ' 2 (l - e 2 ) 1 ' 2 cosz . (28) 

We see that in this special gauge the Delaunay-type equations indeed become canonical, and 
the role of the effective new Hamiltonian is played exactly by the Hamiltonian perturbation 
which emerged earlier in (|14jl. 

Thus we have proven an interesting THEOREM : Even though the gauge-invariant 
Delaunay-type system (J71I - I76|) is not generally canonical, it becomes canonical 
in one special gauge (1221) which we call the "generalised Lagrange gauge." This 
theorem can be proved in a purely Hamiltonian language, as is done in Section 3.3. 

3 Hamilton and Jacobi 

3.1 The concept 

A totally different approach to derivation of the planetary equations is furnished by the tech- 
nique worked out in 1834 - 1835 by Hamilton and refined several years later by Jacobi. In the 
lecture course shaped by 1842 and published as a book in 1866, Jacobi formulated his famous 
theorem and applied it to the celestial motions. Jacobi chose orbital elements that were some 
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combinations of the Keplerian ones. His planetary equations can be easily transformed into 
the Lagrange system by the differentiation chain rule (Subbotin 1968). Later authors used this 
method for a direct derivation of the Lagrange and Delaunay systems, and thus the Hamilton- 
Jacobi approach became a part and parcel of almost any course in celestial mechanics. To some 
of these sources we shall refer below. The full list of pertinent references would be endless, so 
it is easier to single out a couple of books that break the code by offering alternative proofs: 
these exceptions are Kaula (1968) and Brouwer and Clemence (1961). 

Brouwer and Clemence (1961) use the VOP method (like in Lagrange (1808, 1809, 1810)) 
which makes the imposition of the Lagrange constraint explicit. Kaula (1968) undertakes, by 
means of the differentiation chain rule, a direct transition from the Hamilton equations in a 
Cartesian frame to those in terms of orbital elements. As explained in Efroimsky (2002, 2003), 
in Kaula's treatment the Lagrange constraint was imposed tacitly. 

It is far less easy to understand where the implicit gauge fixing is used in the Hamilton- 
Jacobi technique. This subtlety of the Hamilton- Jacobi method is so well camouflaged that 
through the century and a half of the method's life this detail has never been brought to light 
(at least, in the astronomical literature). The necessity of such a constraint is evident: one has 
to choose one out of infinitely many sets of orbital elements describing the physical trajectory. 
Typically, one prefers the set of orbital elements that osculates with the trajectory, so that the 
physical orbit be always tangential to the instantaneous ellipse, in the case of bound orbits, or 
to the instantaneous hyperbola, in the case of flybys. This point is most easily illustrated by the 
following simple example depicted on Fig.l. Consider two coplanar ellipses with one common 
focus. Let both ellipses rotate, in the same direction within their plane, about the shared focus; 
and let us assume that the rotation of one ellipse is faster than that of the other. Now imagine 
that a planet is located at one of the points of these ellipses' intersection, and that the rotation 
of the ellipses is such that the planet is always at the instantaneous point of their intersection. 
One observer will say that the planet is swiftly moving along the slower rotating ellipse, while 
another observer will argue that the planet is slowly moving along the fast-rotating ellipse. 
Both viewpoints are acceptable, because one can divide, in an infinite number of ways, the 
actual motion of the planet into a motion along some ellipse and a simultaneous evolution of 
that ellipse. The Lagrange constraint singles out, of all the multitude of evolving ellipses, 
that unique ellipse which is always tangential to the total (physical) velocity of the body. This 
way of gauge fixing is natural but not necessary. Besides, as we already mentioned, the chosen 
gauge (HJ) will not be preserved in the course of numerical computations. (Sometimes osculating 
elements do not render an intuitive picture of the motion. In such situations other elements 
are preferred. One such example is a circular orbit about an oblate planet. The osculating 
ellipse precesses with the angular velocity of the satellite, and its eccentricity is proportional 
to the oblateness factor J 2 . Under these circumstances the so-called geometric elements are 
more convenient than the osculating ones (Borderies & Longaretti 1987).) 

We remind the reader that the Hamilton- Jacobi treatment is based on the simple facts 
that the same motion can be described by different mutually interconnected canonical sets 
(q, p, TC(q,p) ) and (Q, P, H*(Q,P) ) , and that fulfilment of the Hamilton equations along 
the trajectory makes the infinitesimally small quantities 



d6 = pdq — Hdt 



(29) 



and 



dO = PdQ - H*dt 



(30) 
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perfect differentials. Subtraction of the former from the latter shows that their difference, 

- dW = d6 - dQ = PdQ - pdq - (H* - H) dt (31) 

is a perfect differential, too. Here q, p, Q, and P contain N components each. If we start 
with a system described by (q, p, T~C(q,p) ) , it is worth looking for such a re-parametrisation 
(Q, P, TC*(Q,P)) that the new Hamiltonian H* is constant in time, because in these 
variables the canonical equations simplify. Especially convenient is to find a transformation 
that nullifies the new Hamiltonian TC* , for in this case the new canonical equations will render 




Fig.l. These two coplanar ellipses share one of their foci and are assumed to rotate 
about this common focus in the same direction, always remaining within their plane. 
Suppose that the rotation of one ellipse is faster than that of the other, and that a 
planet is located at one of the points of these ellipses' intersection, P, and that the 
rotation of the ellipses is such that the planet is always at the instantaneous point 
of their intersection. We may say that the planet is swiftly moving along the slower 
rotating ellipse, while it would be equally legitimate to state that the planet is slowly 
moving along the fast-rotating ellipse. Both interpretations are valid, because one 
can divide, in an infinite number of ways, the actual motion of the planet into a 
motion along some ellipse and a simultaneous evolution of that ellipse. The Lagrange 
constraint (j3J) singles out, of all the multitude of evolving ellipses, that unique ellipse 
which is always tangential to the total, physical, velocity of the planet. 
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the variables (Q, P) constant. One way of seeking such transformations is to consider W 
as a function of only q , Q , and t . Under this assertion, the above equation will entail: 

dW dW dW 

— -5— dt — —— dQ — — — dq = P dQ - pdq + (H - H*) dt (32) 

at oQ dq 

whence 

dW dW dW 

P = SQ ' »=Hi ' H + lH= H ' ' < 33 > 
The function W can be then found by solving the Jacobi equation 

/ dW \ dW 

"('•W t ) + -sr = Tr (34) 

where 7i* is a constant. It is very convenient to make it equal to zero. Then the above 
equation can be easily solved in the unperturbed (reduced) two-body setting. This solution, 
which has long been known, is presented, in a very compact form, in Appendix 2. It turns 
out that, if the spherical coordinates and their conjugate momenta are taken as a starting 
point, then the eventual canonical variables Q, P corresponding to H*(Q, P) = are the 
Delaunay elements: 

Qi = L — y/JTa , Pi = - M 



Q 2 = G = yjna{\- e 2 ) , P 2 = - to , (35) 

Q 3 = H = \j[ia (1 — e 2 ) cos i , P 3 = — 
3.2 Where can free cheese be found? 

The transition from two-body to N-body celestial mechanics is presented in numerous books. 
However, none of them explain how the Lagrange constraint is implicitly involved in the for- 
malism. 

Before we move on, let us cast a look back at what has been accomplished in the two-body 
case. We started out with a Hamiltonian problem ( q, p, Ti. ) and re-formulated its equations 
of motion 

. on . m 

q = ' P = ~ (36) 

in terms of another set (Q, P, 7i* ) : 

q = 0(Q, P, t) 

(37) 

p = P, t) 

so that the above equations are mathematically equivalent to the new system 

Q ~ ~W ' p ~ " w ■ (38) 

The simple nature of the 2-body setting enabled us to carry out this transition so that our new 
Hamiltonian Ti* vanishes and the variables Q and P are, therefore, constants. This was 
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achieved by means of a transformation-generating function W(q, Q, t) obeying the Jacobi 
equation (|34|). After formula ()89|) for this function was written down, the explicit form of 
dependence (|57)l can be found through the relations P = — dW/ dQ . This is given by (p2l . 

To make this machinery function in an N-body setting, let us first consider a disturbed 
two-body case. The number of degrees of freedom is still the same (three coordinates q and 
three conjugate momenta p), but the initial Hamiltonian is perturbed: 

. = din + ah) = _ din + ah) 

dp ' dq 

Trying to implement the Hamilton- Jacobi method (}3*2*j) - ([34)1 . for the new Hamiltonians (H + 
AH) , (H* + AH) and for some generating function V(q, Q, t) , we shall arrive at 

dV dV dV 
~^dt - — dQ - — dq = PdQ - pdq + [iH + AH) - (H* + AH) ] dt , 40 
dt dQ dq 

dV dV dV 

p = -1>q • P = W ' U + AH+ m =H " + AH - (41) 

( dV \ dV 

+ 1* =n " ■ < 42 > 

We see that V obeys the same equation as W and, therefore, may be chosen to coincide 
with it. Hence, the new, perturbed, solution (q, p) will be expressed through the perturbed 
"constants" Q(t) and P(t) in the same fashion as the old, undisturbed, q and p were 
expressed through the old constants Q and P : 

q = <f)(Q(t), P(t), t) , 

(43) 

p = ^{Qit), P(t), t) , 

<f) and ip being the same functions as those in (JHTj) . Benefitting from this form-invariance, one 
can now master the N-particle problem. To this end, one should choose the transformation- 
generating function V to be additive over the particles, whereafter the content of subsection 
3.1 shall be repeated verbatim for each of the bodies, separately. In the end of this endeavour 
one will arrive to iV — 1 Delaunay sets similar to except that now these sets will be 

constituted by instantaneous orbital elements. The extension of the preceding subsection's 
content to the N-body case seems to be most straightforward and to involve no additional 
assumptions. To dispel this illusion, two things should be emphasised. One, self-evident, fact 
is that the quantities Q and P are no longer conserved after the disturbance AH is added 
to the zero Hamiltonian H* . The second circumstance is that a change in a Hamiltonian 
implies an appropriate alteration of the Lagrangian. In the simple case of AH being a 
function of the coordinates and time only (not of velocities or momenta), addition of AH to 
the Hamiltonian implies addition of its opposite to the Lagrangian. Since this extra term bears 
no dependence upon velocities, the expressions for momenta through the coordinates and time 
will stay form-invariant. Hence (if the Lagrangian is not singular), the functional dependence 
of the velocities upon the coordinates and momenta will, also, preserve their functional form 
v(q, p, t) : 

• t 7 9C(q, q,t) 

without perturbation: p = q — v(q,p,t) , 

dq 



12 



with perturbation : p = - - —Ql — — ^— ^ — — — = — — - =>- q = v(q, p, t) . 

oq oq 

where the new, perturbed, dependence q = v[q(Q(t), P(t), t), p(Q(t), P(t), t), t] has the 
same functional form as the old one, q = v[q(Q, P, t), p(Q, P, t), t] . Together with (jJ3J), 
this means that the dependence of the new q upon the new P(t) and Q(t) will have the 
same functional form as the dependence of the old q upon the constants Q and P : 

j t q(Q(t),P(t),t) = | q(Q(t),P(t),t) . (45) 



In other words, 



£gj|-A = 0, (46) 



=1 



where Di denotes the set of perturbed variables (Q(t), P(t)) . In the astronomical applica- 
tions, Di may stand for the Delaunay set. 

This is the implicit condition under which the Hamilton- Jacobi method works (in the above 
case of velocity-independent disturbance). Violation of (J46|) would invalidate our cornerstone 
assumption (|38|). This circumstance imposes a severe restriction on the applicability of the 
Hamilton- Jacobi theory. In the astronomical context, this means that the Delaunay elements 
(|92*|) must be osculating. Indeed, if Di denote a set of orbital elements, then expression (j4*Uj) 
is equivalent to the Lagrange constraint (JH) discussed in Section 1. There the constraint was 
imposed upon the Keplerian elements; however, its equivalence to (|4fij) . which is written in 
terms of the Delaunay variables, can be easily proven by the differentiation chain rule. 



(44) 



3.3 The case of momentum-dependent disturbances 

When the perturbation of the Lagrangian depends also upon velocities (and, therefore, the 
Hamiltonian perturbation carries dependence upon the canonical momenta), the special gauge 
(j22J) wherein the Delaunay- type system preserves its canonicity differs from the Lagrange gauge. 
This was proven in subsection 2.3 in the Lagrangian language. Now we shall study this in 
Hamiltonian terms. Our explanation will be sufficiently general and will surpass the celestial- 
mechanics setting. For this reason we shall use notations q, p , not r, p . The development 
will, as ever, begin with an unperturbed system described by canonical variables obeying 

. = m . = _ m 

^ dp ' ^ dq 
This dynamics may be re-formulated in terms of the new quantities (Q, P) : 

q = <KQ, P, t) 

(48) 

p = P, t) 

so that the Hamiltonian equations ()47|) are equivalent to 

q = P = - d 2L (49) 

W dP ' dQ ■ 1 ' 
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For simplicity, we shall assume that H* is zero. Then the new canonical variables will play 
the role of adjustable constants upon which the solution P5|l of (jj7j) depends. 
We now wish to know under what circumstances a modified canonical system 

d(H + AH) . d(H + AH) 
q = ^ " , P = ~ -^—g-q > An = An (^ P. *) ( 5 °) 



will be satisfied by the solution 



q = cf>{Q{t), P{t), t) 



(51) 



p = iP(Q(t), P(t), t) 
of the same functional form as PHj) but with time-dependent parameters obeying 



• _ dAH p _ _d_AH 

Q ~ 8P ' dQ • [b } 

This situation is of a more general sort than that addressed in subsection 3.2, in that the 
perturbation AH now depends also upon the momentum. 

Under the assumption of being (at least, locally) invertible, substitution of the equal- 
ities 

• dAH dAH dq_ dAH dp_ 

Q ~ 8P ~ dq dP + dp dP ^ ' 



and 



dQ ~ dq~ dQ dp~ dQ ^ ' 



into the expression for velocity 



dq dq_ d dq_ p 

q ~ dt + dQ Q + dP F (55j 



leads to 



dq I dq dq dq dq\ dAH ( dq dp dq dp\ dAH , , 

q = dt + [dQ dP ~ dPdQJ dq + [dQ dP ~~ dP dQ J dp ' ( } 

Here the coefficient accompanying dAH/dq identically vanishes, while that accompanying 
dAH I dp coincides with the Jacobian of the canonical transformation and is, therefore, unity: 

dq_ dp_ dq_dp_ 

dQ dP dP dQ 1 ' 

So if we introduce, in the spirit of (JBJ), notation 



then (jSUjl will lead to 



(dAH\ 

3 + M, ■ (59) 



q,t 
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Expression ([59)1 establishes the link between the regular VOP method and the canonical treat- 
ment. It shows that, to preserve the symplectic description, one must always choose a particular 
gauge $ = d ATi/dp . Needless to say, this is exactly the generalised Lagrange gauge (J2*2*|l 
discussed in subsection 2.3. A direct, though very short, proof is as follows. 

On the one hand, the Hamilton equation for the perturbed Hamiltonian (|14|) is: 

. d (H + AH) d AH 

q = ~ = p + -t, — (60) 

op op 

while, on the other hand, the definition of momentum entails, for the Lagrangian ([120: 

d(C(q, q, t) + AC(q, q, t)) . d AC 
P = dq = q + ~W • (61) 

By comparing the latter with the former we arrive at: 

* - PTL - - (^L 

which coincides with (12 2|) . Thus we see that transformation (J48|) being canonical is equivalent to 
the partition of the physical velocity q in a manner prescribed by ()59j) . where $ = d ATi/dp . 
This is equivalent to our Theorem from subsection 2.3. Evidently, for disturbances dependent 
solely upon the coordinates, we return to the case explained in subsection 3.2 (equations (f4"Sl 
- BSj)): in that case, the Hamiltonian formulation of the problem demanded imposition of the 
Lagrange constraint (|46j) . 

To draw to a close, the generalised Lagrange constraint, $ = — dAC/dq , is stiffly 
embedded in the Hamilton- Jacobi technique. Hence this technique breaks the gauge invariance 
and is unfit (at least, in its straightforward form) to describe the gauge symmetry of the 
planetary equations. It is necessary to sacrifice gauge freedom by imposing the generalised 
Lagrange constraint to make use of the Hamilton- Jacobi development. 

In this special gauge, the perturbed momentum coincides with the unperturbed one (which 
was equal to g ). Indeed, we can rewrite (JBTj) as 

_ d(C(q, q, t) + AC(g, q t) ) . 
p = — = q - = g (63) 

which means that, in the astronomical implementation of this theory, the Hamilton- Jacobi 
treatment necessarily demands the orbital elements to osculate in the phase space. Naturally, 
this demand reduces to that of regular osculation in the simple case of velocity-independent 
AC that was explored in subsection 3.2. 



4 Conclusions 

We have studied, in an arbitrary gauge, the VOP method in celestial mechanics in the case 
when the perturbation depends on both positions and velocities. Such situations emerge when 
relativistic corrections to the Newton law are taken into account or when the VOP method 
is employed in noninertial frames of reference (a satellite orbiting a precessing planet being 
one such example). The gauge-invariant (and generalised to the case of velocity-dependent 
disturbances) Delaunay-type system of equations is not canonical. We, though, have proven 
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a theorem establishing a particular gauge (which coincides with the Lagrange gauge in the 
absence of velocity dependence of the perturbation) that renders this system canonical. We 
called that gauge the "generalised Lagrange gauge." 

We have explained where the Lagrange constraint tacitly enters the Hamilton- Jacobi deriva- 
tion of the Delaunay equations. This constraint turns out to be an inseparable (though not 
easily visible) part of the method: in the case of momentum-independent disturbances, the 
N-body generalisation of the two-body Hamilton- Jacobi technique is legitimate only if we use 
orbital elements that are osculating, i.e., if we exploit only the instantaneous ellipses (or hy- 
perbolae, in the flyby case) that are always tangential to the velocity vector. Oddly enough, 
an explicit mention of this circumstance has not appeared in the astronomical literature (at 
least, to the best of our knowledge). 

In the case of momentum-dependent disturbances, the above restriction generalises, in 
that the instantaneous ellipses (hyperbolae) must be osculating in the phase space. This is 
equivalent to the imposition of the generalised Lagrange gauge. 

Comparing the good old VOP method with that based on the Jacobi theorem, we have 
to acknowledge that the elegance of the latter does not outweigh the power of the former. If 
we decide to explore the infinite multitude of gauges or to study the numerical-error-invoked 
gauge drift, we shall not be able to employ the Hamilton- Jacobi theory without additional 
structure. However, the direct VOP method unencumbered with the canonicity demand will 
immediately yield gauge-invariant equations for the Delaunay elements obeying an arbitrary 
gauge condition 

^ df dDj * . . . . 

^mnr = *> • (64 > 

$ being some function of time and elements Di . In Efroimsky (2002) these equations were 
written down for the case of velocity-independent perturbation. If the disturbing force depends 
also upon velocities, the Delaunay-type equations will acquire even more terms and will read 
as (J7T1 - I7BJ) . In the simple case of a velocity-independent disturbance, any supplementary con- 
dition different from that of Lagrange will drive the Delaunay system away from its canonical 
form. If we permit the disturbing force to depend also upon velocities, the Delaunay equations 
will retain their canonicity only in the generalised Lagrange gauge. 

In the language of modern physics this may be put in the following wording. N-body 
celestial mechanics is a gauge theory but is not genuinely symplectic insofar as the language 
of orbital elements is used. It, though, becomes canonical in the generalised Lagrange gauge. 

The applications of this formalism to motions in non-inertial frames of reference will be 
studied in Efroimsky & Goldreich (2003). Some other applications were addressed in Slabinski 
(2003). 



Acknowledgments 

The authors are grateful to Jean Kovalevsky for an e-mail conversation, which he had with 
one of us (ME) and which initiated this work. ME would like to acknowledge his numerous 
stimulating discussions with Marc Murison, William Newman, and Victor Slabinski. Research 
by ME was supported by NASA grant W-19948. Research by PG was partially supported by 
NSF grant AST 00-98301. 



16 



Appendix 1: 



Gauge- invariant equations of Lagrange and Delaunay types 

We present the gauge-invariant Lagrange-type equations. They follow from (JT§j) if we take into 
account the gauge- invariance of matrix [CjCj] defined by ©• We denote by AH the pertur- 
bation of the Hamiltonian, connected through (jl4j) with that of the Lagrangian. The latter, 
in its turn, is connected through (fT7)|) with the disturbing force (and acts as the customary 
disturbing function when the perturbations are devoid of velocity dependence). 



da 
~dt 



2 

na 



d(-AH) dAC d 



dM 



$ + 



dr 0M o 



dAC\ d$ 



3> 



3M' 
dr , 



8i d 



Or J dM dM dt 



$ + 



dAC" 
dr , 



(65) 



de 
dt 



1-e 2 



na* e 



d ( - AH) dAC d_ f - dAC " 
dM 5f da 1 d'r , 



$ + 



dAC\ dt 



dr ) dM Q dM dt 



df d ( z dAC" 
dr , 



(66) 



,2U/2 



na? e 



d(-AH) dAC d dAC" 



duo 



dr du 



dr 



dAC\ dg df d ( =? dAC\ 



dr J du duj dt 



dr J 



du 



— cos i 



dt na 2 (l — e 2 ) 1 / 2 sin i 



d ( - AH) _ dAC d_ ( - dAC " 
di d'r di\ d'r , 



$ + 



dAC \ dg _ df_ d_ / - dAC " 
<9F / di di dt\ d'r , 



+ 



(67) 



1-e 2 ) 1 / 2 



na? e 



d ( - AH) 

de 



dAC d (* dAC" 
$ + 



dr de 



dr 
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d? / de de dt \ d'r ) 



di cos i 

dt na 2 {l — e 2 ) 1 / 2 sin i 



d(-AH) dAC d d AC 



duo 



dr duo 



dr 



$ + 



dAC\ dg di d /-j dAC\ 



dr 



y ) duo duo dt 



dr 



V) 



na 2 (l — e 2 ) 1 / 2 sin i 



d(-AH) dAC d dAC 



dtt 



dr dtt V dr 



(68) 



' - dA£\ dg df d_ U dAC \ 



dVL 1 

dt na 2 (l — e 2 ) 1 / 2 sin i 



d(-AH) dAC d dAC 



di 



dr di 



dr 



$ + 



dAC\ dg di d dAC\ 



dr 



? ) di di dt [® + d'r ) 



(69) 



dM Q 
~df 



1 - e 2 
no? e 



d(-AH) dACd_(^ + dAC 



de 



dr de 



dr 



$ + 



<9A£\ <9g 9f d <9A£ 



<9r 



r / de de dt 



$ + 



<9r 



(70) 



2 

na 



d(-AH) dAC d (~ , <9A£ 



<9a 



<9r 



' - + 9A£\ 9g _ 9f d f % + d_AC 
K dr J da da dt \ dr 
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Similarly, the gauge-invariant Delaunay-type system can be written down as: 

dL _ d{ — AH) dAL d d AC\ d AC\ dg df d ( - d AC\ 

~dt~ dM dfdMo\ + ~dYj~[ + ~dF) dM~o ~ dW Jt \ + ~W ) 

dM, = d(-AH) dAC d f $ dAC\ ( $ dAC\ d£ flf d_ ( $ d AC \ 

dt dL df dL \ d? ) \ d'v ) dL dL dt\ d? ) 



dG d ( - AH) _ dAC d_ / - d AC \ ( - dA£\ gg df d_ / - d AC ^ 

dt duj dr du \ d 'v J \ <9r / <9u; <9a> ^ \ <9r y 

rf^ g ( - A7i) <9A£ g dA£\ dAC \ dg df_ d_ ( - 3A£ X 

d* 7 ~ + ~dFdG { + d'v ) + V + a?/5G + 5G^V + dF , 

dff ^ ( - A7i) <9A£ g <9A£ \ dA£ \ dg df <9A£\ 




rffi = d(-AH) dAC d dAC\ dAC\ dg A* 

d# d? dH\ dt J \ d'v ) dH dH dt 1 

where 

L = ^a 1 ' 2 , G = ^ 2 a 1 ' 2 (l - e 2 ) 1 ' 2 , H = ^ 2 a 1 ' 2 (l - e 2 ) 1 ' 2 cos i (77) 

— * — * _ 

and the symbols f, g denote the functional dependencies of the gauge, position and velocity 
upon the Delaunay, not Keplerian elements, and therefore these are functions different from 
<fr, f, g used in (j65l - 170|) where they stood for the dependencies upon the Kepler elements. 
(In Efroimsky (2002) the dependencies f, g upon the Delaunay variables were equipped 
with tilde, to distinguish them from the dependencies upon the Kepler coordinates.) 

The above equations do not merely repeat those derived earlier in Efroimsky (2002, 2003), 
but generalise them to the case of a perturbation AC which is both position- and velocity- 
dependent. For this reason, our gauge-invariant equations can be employed not only in an 
inertial frame but also in a wobbling one. 

To employ the gauge-invariant equations in analytical calculations is a delicate task: one 

— *■ 

should always keep in mind that, in case $ is chosen to depend not only upon time but also 
upon the "constants" (but not upon their derivatives), the right-hand sides of these equation 
will implicitly contain the first derivatives dCi/dt , and one will have to move them to the 
left-hand sides (much like in the transition from (|10|) to 

Appendix 2: 

The Hamilton-Jacobi method in celestial mechanics 



The Jacobi equation (J53jl is a PDE of the first order, in (N+l) variables (q n , t) , and its 
complete integral W(q, Q, t) will depend upon N+l constants a n (Jeffreys and Jeffreys 
1972, Courant and Hilbert 1989). One of these constants, a N+ i , will be additive, because 
W enters the above equation only through its derivatives. Since both Hamiltonians are, too, 
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defined up to some constant / , then the solution to ()34|) must contain that constant multiplied 
by the time: 



W(q, ai, ... , a N , a N+1 , t) = W(q, a x , ... , a N , t) - (t - t Q ) /(ai, ... , a N ) 

(78) 

= W(q, ax, ... , a N , t) - t f(a ly ... , a N ) - a N+1 

where the fiducial epoch is connected to the constants through t a = —a N+ \/f , and the 
function W depends upon N constants only. As the total number of independent adjustable 
parameters is N + 1 , the constant / cannot be independent but must rather be a function 
of ai, ... , a N , a N+ i . Since we agreed that the constant a N+ i is additive and shows itself 
nowhere else, it will be sufficient to consider / as a function of the rest N parameters only. 
(In principle, it is technically possible to involve the constant a N+ i , i.e., the reference epoch, 
into the mutual transformations between the other constants. However, in the applications 
that we shall consider, we shall encounter only equations autonomous in time, and so there 
will be no need to treat a N+ i as a parameter to vary. Hence, in what follows we shall simply 
ignore its existence.) The new function W obeys the simplified Jacobi equation 

/ dW(q,a 1 ,...,a N ,t) \ dW(q, m, ... , a N , t) 
H \q , — , 1 1 + — = f(ai, ... , a N ) + H (79) 

As agreed above, 7i* is a constant. Hence, we can state about this constant all the same as 
about the constant / : since the integral W can contain no more than iV + 1 adjustable 
parameters a 1; ... , a N , a N+ i , and since we ignore the existence of a N+ i , the constant 7i* 
must be a function of the remaining N parameters: TC* = 7i*{a\, ... , a N ). 

Now, in case TC depends only upon (q, p) and lacks an explicit time dependence, then 
so will W; and the above equation will very considerably simplify: 

/ dW(q, ai, a N )\ . 
n\q, 7^ I = f{a u ... , a N ) + H (a u ... , a N ) , (80) 

where we deliberately avoided absorbing the constant Hamiltonian 7i* into the function / . 

Whenever the integral W can be found explicitly, the constants {a\, ... , a N ) can be 
identified with the new coordinates Q , whereafter the new momenta will be calculated through 
P = — dWj OQ . In the special case of zero 7i* , the new momenta become constants, because 
they obey the canonical equations with a vanishing Hamiltonian. In the case where 7i* is a 
nonzero constant, it must, as explained above, be a function of all or some of the independent 
parameters (ai, ... , a N ) , and therefore, all or some of the new momenta P will be evolving 
in time. 

Since it is sufficient to find only one solution to the Jacobi equation, one can seek it by 
means of the variable-separation method: equation (|80p will solve in the special case when the 
generating function (f7%|) is separable: 

N 

W(q!, ... , q N , ai, ... , a N ) = W * > Q i> - > a Jv) ( 81 ) 

This theory works very well in application to the unperturbed (two-body) problem of 
celestial mechanics, a problem that is simple due to its mathematical equivalence to the grav- 
itationally bound motion of a reduced mass rn p i anet m sun / (m p i anet + m sun ) about a fixed 
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centre of mass vn v \ anet + m sun . If one begins with the (reduced) two-body Hamiltonian in 
the spherical coordinates 

qi = r , q 2 = (f> , q 3 = 9 (82) 
(where x = r cos cos 6 , y = r cos <fi sin 9 , z = r sin 0) , then the expression for Lagrangian, 

L = T - II = \ { qi f + \ { qi f (q 2 f + \ { qi f {q 3 f cos 2 g 2 + ^ (83) 
L L L qi 

will yield the following formulae for the momenta: 

Pi = = 9i , P2 = W7- = q x q2 , Ps = W7- = q x q 3 cos g 2 (84) 
oqi oq 2 oq 3 

whence the initial Hamiltonian will read: 



n = J2pi - £ = \vl + ikvl + 9 2 1 2 ^3 - - • (85) 

2 2gf 2gf cos 2 g 2 <Zi 



Then the Hamilton- Jacobi equation (30) will look like this: 

1 (dW\ 2 1 (dW\ 2 1 (dW\ 2 /I dW 

2 \~dq~i) + 23 \dq~2) + 2g 2 cos 2 g 2 I^J ~ J x ~ ~dt ~ ' (86) 
while the auxiliary function H 7 defined through (fTSj) will obey 

l(fV + 1 1 (f) 2 _ii _/_„.-„. (87, 

2 V 9gi / 2gf \ <9g 2 / 2gf cos 2 q 2 \ dq 3 J q 1 

A lengthy but elementary calculation (presented, with some inessential variations, in Plummer 
1918, Smart 1953, Pollard 1966, Kovalevsky 1967, Stiefel and Scheifele 1971, and many other 
books) shows that, for a constant H* and in the ansatz (j8Tf . the integral of (}80|) takes the 
form: 

W = Wi (qi, oi, a 2 , a 3 ) + W 2 (q 2 , a t , a 2 , a 3 ) + W 3 (q 3 , a x , a 2 , a 3 ) = 

(88) 

9lit) e x (2(f + 7C) + ^ - 4V ^qi + f* e 2 (a\ - ^f-V \q 2 + Ta 3 dq 3 

9i (<<,) V qi qi) Jo \ cos 2 g 2 y jo 

where the epoch and factors ei, e 2 may be taken as in Kovalevsky (1967): time t Q is the instant 
of perigee passage; factor e\ is chosen to be + 1 when q x = r is increasing, and is — 1 when r 
is decreasing; factor e 2 is + 1 when g 2 = is increasing, and is — 1 otherwise. This way the 
quantities under the first and second integration signs have continuous derivatives. To draw 
conclusions, in the two-body case we have a transformation-generating function 

W = W + tf(a u ...,a N ) = 

(89) 

9l(t) £i (2(f + H*) + — — 4V "dq, + f* e 2 [a\ - * ' dq 2 + f a 3 dq 3 + tf 
91 (to) V Qi QiJ Jo \ cos z q 2 J Jo 
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whose time-independent component W enters equation (|80p. The first integration in (|89|) 
contains the functions f(a\, ... , a N ) and H*(ai, ... , a N ) , so that in the end of the day W 
depends on the N constants ai, ... , a N (not to mention the neglected t D , i.e., the a N+ i ). 

Different authors deal differently with the sum (/ + TC*) emerging in (|89|) . Smart (1953) 
and Kovalevsky (1967) prefer to put 

/ = , H* = a x , a 1 = -fi/{2a) , (90) 

whereupon the new momentum Pi — — dW/ 8Q\ = — dW/ da\ becomes time-dependent 
(and turns out to equal —t + t a .) An alternative choice, which, in our opinion, better reflects 
the advantages of the Hamilton- Jacobi theory, is furnished by Plummer (1918): 

/ = a,\ , 7i* = , a% — y/JTa . (91) 

This entails the following correspondence between the new canonical variables (the Delaunay 
elements) and the Keplerian orbital coordinates: 

Qx = ai = y/JTa , Pi = - M Q , 

Q 2 = a 2 = y/fj,a(l- e 2 ) , P 2 = - uj , (92) 



Q 3 = a 3 = sj fia (1 - e 2 ) cos i , P 3 = - Q 

Everywhere in this paper we follow the convention (|91j) and denote the above variables Qi, Q 2 , Q3 
by L, G, H, correspondingly (as is normally done in the astronomical literature) 
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